Population imbalanced fermions in harmonically trapped optical lattices 
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The attractive Fermi-Hubbard Hamiltonian is solved via the Bogoliubov-de Gennes formalism to 
analyze the ground state phases of population imbalanced fermion mixtures in harmonically trapped 
two-dimensional optical lattices. In the low density limit the superfluid order parameter modulates 
in the radial direction towards the trap edges to accommodate the unpaired fermions that are pushed 
away from the trap center with a single peak in their density. However in the high density limit while 
the order parameter modulates in the radial direction towards the trap center for low imbalance, 
it also modulates towards the trap edges with increasing imbalance until the superfluid to normal 
phase transition occurs beyond a critical imbalance. This leads to a single peak in the density of 
unpaired fermions for low and high imbalance but leads to double peaks for intermediate imbalance. 
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The phase diagram of dilute population imbalanced 
fermion mixtures has been recently studied showing su- 
perfluid to normal phase transition with increasing im- 
balance as well as a phase separation between paired and 
unpaired fermions [H, i> i, El]. These recent works are 
extentions of the earlier works on dilute population bal- 
anced mixtures where a crossover from Bardeen-Cooper- 
Schrieffer (BCS) to Bose-Einstein condensation (BEC) 
type superfluidity is observed as a function of the attrac- 
tive fermion-fermion interaction strength [H, H H i, i] ■ 

Arguably understanding the phase diagram of fermion 
mixtures in optical lattices is one of the next frontiers 
in cold atoms research because of their great tunability. 
In addition to the particle populations and the particle- 
particle interaction strengths, one can also precisely con- 
trol the particle tunnelings, the lattice dimensionality 
and the lattice geometry. For instance experimental evi- 
dence for the superfluid and the insulating phases of pop- 
ulation balanced mixtures have been recently reported in 
trapped optical lattice s II CI. after overcoming some ear- 
lier difficulties [ll|, EH Il3lll4| . This recent work has also 
opened the possibility of studying many-body properties 
of population imbalanced mixtures in optical lattices. 

Earlier theoretical works on population imbalanced 
fermion mixtures i n op tical lattices were limited to ho- 
mogenous systems [H[ HU , and they showed rich phase 
diagrams involving BCS type nonmodulating and Fuldc- 
Ferrell-Larkin-Ovchinnikov (FFLO) type spatially mod- 
ulating superfluid phases in addition to insulating and 
normal phases. Furthermore the phase diagram of popu- 
lation imbalanced mixtures in harmonically trapped op- 
tical lattices has been recently discussed within the semi- 
classical local density approximation (LDA) [l7j]. How- 
ever it has been previously shown that the LDA type 
methods are not sufficient to describe even the dilute 
population imbalanced mixtures without an optical lat- 
tice [HI, [HI . In this manuscript we therefore analyze the 
ground state phases of fermion mixtures in harmonically 
trapped two-dimensional optical lattices via using the 
fully quantum mechanical Bogoliubov-de Gennes (BdG) 



method where the trapping potential is included exactly 
at the mean-field level. 

Our main results are as follows. In the low density 
limit the superfluid order parameter modulates in the 
radial direction towards the trap edges to accommodate 
the unpaired fermions that are pushed away from the 
trap center with a single peak in their density. These 
fin ding s are in good agreement with the recent theoreti- 
cal [IlOjIlIl and experimental [2, Hi, El 

findings on di- 
lute population imbalanced mixtures without an optical 
lattice. However in the high density limit while the order 
parameter modulates in the radial direction towards the 
trap center for low imbalance, it also modulates towards 
the trap edges with increasing imbalance until the super- 
fluid to normal phase transition occurs beyond a critical 
imbalance. This leads to a single peak in the density of 
unpaired fermions for low and high imbalance but leads 
to double peaks for intermediate imbalance. 

BdG formalism: To achieve these results we solve the 
Fermi-Hubbard Hamiltonian 



Hfh = - ^ ' <'*\, r <." j ", r .T"<- 
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where a\ a (a iiCr ) creates (annihilates) a pseudo-spin a 
fermion at lattice site i, tij, a and Uij > are the 
particle-particle tunneling and the density-density inter- 
action matrix elements, \i a is the chemical potential, and 
Vi.a = a CT |ri| 2 /2 is the trapping potential at position 
with a a = rn a u>0 such that the trapping potential is cen- 
tered at the origin. Here the label a identifies T or J, 
fermions and allows a fermions to have equal or unequal 
masses controlled by ii.j l(T and/or to have equal or un- 
equal populations controlled by \i a . 

In the mean-field approximation for the superfluid 
phase, the Fermi-Hubbard Hamiltonian reduces to 
H = -T, i ,j,a t hj^4,a a j^ ~ E,vO* - V ita )al a a ita - 



t t + A^a^TOjj - [At^/Uij), where the 
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self-consistent field Ajj = Uij(ai^aj_i) is the su- 
pcrfluid order parameter and (...) is a thermal aver- 
age. The mean-field Hamiltonian can be diagonal- 
ized via the Bogoliubov-Valatin transformation ai^ — 
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T ], where 7* (-y n>(T ) cre- 



ates (annihilates) a pseudo-spin cr quasiparticle with the 
wavefunction u n _i^ a (v n ^ -a ), and ,S| = +1 and = — f. 
This leads to the BdG equations 



E 



T; 



A* 
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(2) 



where T it j a = —tij tC — (/i^ — Vi !<T )di t j is the diagonal el- 
ement and <5i.j is the Kronecker delta. Here e„ j(T > 
are the eigenvalues and (p n ,i.<7 are the eigenfunctions 
given by (f a ,i,^ = iK,i,P for the T and fa,i,^ = 

{Vn,i,h ~ u n,i,i) for the J. eigenvalues. Since solutions to 
the BdG equations are invariant under the transforma- 
tion Vn ti< i -> < i4iT , Wn.i.J. -> ~<,a and e n,l ~ e «,T' 
it is sufficient to solve only for u n ^ = u n ^, v n ^ = v Ui i t i 
and e n = t n ^ as long as we keep all the solutions with 
positive and negative e„. 

In Eq. © the superfluid order parameter Ajj is 
given by A iti = ~J2n U i,j u nAKjf( e n) where f(x) = 
f/[exp(x/T) + 1] is the Fermi function and T is the 
temperature. Notice that this equation is free from the 
ultraviolet divergence and therefore it does not explic- 
itly involve any energy cut-off since the lattice spac- 
ing provides an implicit cut-off. Eq. ^ and the order 
parameter equation have to be solved self-consistently 
with the number equations < ni y(T < 1 = (a\ a a.i t a) 
for the a fermions, where = \u nj i\ 2 f(e n ) and 
n U = J2 n \ v n,i\ 2 f(- e n) such that N a = J2i n h°- deter- 
mines /i CT . In the following we consider only attractive 



and onsite s-wave interactions and set L\ 
Un > 0. Notice that this leads to Ai 



UoSij with 
Ai<5i.j. Further- 



more fermions are allowed to tunnel only to the nearest 
neighbor sites and thus t^j^ — t a 5ij±\. 

Ground state phases: We now analyze the ground state 
phases of fcrmion mixtures with equal masses (too — 
?7ij = mjj, equal tunnclings (to = t-\ = i|) and equal 
trapping potentials (ao = a-\ = a\) but with unequal 
chemical potentials. The theoretical parameters to and 
Uo can be expressed in terms of the experimental pa- 
rameters of the two-dimensional optical lattice poten- 
tial VL(x,y) = VL[sin 2 (irx/a) + sin 2 (-7r?//a)] via the re- 
lations [2| t = (4^ r /V^)(VL/S r ) 3 / 4 cxp(-2 v /TV^) 
and Uo = — V&KaFE r (VL/E r ) 3 / 4 /a. Here a is half of 
the laser wavelength which corresponds to the lattice 
spacing, Vl is the depth of the optical lattice potential, 
E r = h 2 Tr 2 /(2tooo 2 ) is the recoil energy and of is the 
two-body scattering length in vacuum. The experimen- 
tal parameters Vl, a and aF can be tuned by varying the 
laser intensity, the laser wavelength and the externally 
applied magnetic field via using the Feshbach resonances, 
respectively, which makes optical lattices ideal systems to 
simulate the Fermi-Hubbard Hamiltonian. 



For numerical purposes the superfluid order parame- 
ter is assumed to be real (A^ — A*). This is sufficient 
to describe the nonmodulating and the spatially modu- 
lating superfluid phases in addition to the normal and 
the band insulator phases. We also take Uo = 3io and 
Vo = aoo 2 /2 = 0.02io as the strength of the weak onsite 
interactions and the weak trapping potentials, respec- 
tively, and perform calculations on a two-dimensional 
square lattice with a length of L = 50a in both directions. 
The trap center is located at r c = (x = 0a,y = Qa). We 
want to emphasize that similar calculations can be also 
performed for three-dimensional optical lattices. How- 
ever they are computationally much more demanding and 
we do not expect any qualitative difference between our 
two-dimensional results and the three-dimensional ones. 

We fix the total number of fermions N = + iVj 
to N ss 270 (corresponding to ^ « Oto) in the low den- 
sity and to N « 1570 (corresponding to /i ps 5to) in the 
high density case where (i = (/if + /ij)/2, while we vary 
the population imbalance P — (TVf — Ni)/N or equiv- 
alently 5fi ~ (/Xj — ^j)/2. For these parameters it is 
important to notice that the trapping potential provides 
a soft boundary leading to a finite system, and there- 
fore it simplifies the numerical calculations considerably 
in comparison to infinite systems. Next we present self- 
consistent solutions of Eq. ([2]) with the order parameter 
and the number equations. 
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FIG. 1: (Color online) We show (a) the order parameter Aj (in 
units of to) and (b) the population difference pi — m,t — n i.l 
(per lattice site) for the low density case as a function of 
distance x (in units of a) from the trap center. Here y = 0a. 

(I) Low density mixtures: In Fig. [1] we show the super- 
fluid order parameter A^ and the population difference 
per lattice site pi — n^-j- — n^i for the low density case 
where N « 270. When U = and 7V T = N L = N/2, 
the maximum filling of this case corresponds to an al- 
most half-filled band with n.; j(T 0.5 at the trap center. 
For such low densities we expect that our results for the 
trapped mixtures with an optical lattice to recover the 
previously obtained results for the tra ppe d dilute mix- 
tures without an optical lattice [H, [H, |20|| . This occurs 
when the interparticle separation becomes much longer 
than a such that the particles do not feel the presence 
of a lattice potential. However to understand the ground 
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state phases of population imbalanced mixtures, it is very 
illustrative to first discuss the population balanced case. 

For a weakly attracting population balanced mixture 
with Uq = 3in and 5/1 = 0, the order parameter A, is 
finite around the trap center for distances |ri| < 15a, 
and therefore the ground state corresponds to a BCS 
type superfluid. For longer distances |ri| > 16a away 
from the trap center, gradually decreases until it 
eventually vanishes when the densities become very low 
Tii t -\ = n^i w 0. These features can be seen in Fig. QJa), 
and they are in good agreement with the earlier exper- 
iments involving population balanced mixtures without 
an optical lattice [5|, la, LD, la, l2( ■ 

In the case of population imbalanced mixtures, we find 
that Aj modulates in the radial direction towards the 
trap edges to accommodate the unpaired fermions. How- 
ever Aj decreases with increasing population imbalance 
as shown in Fig. HJa), and it vanishes entirely beyond a 
critical imbalance signaling a transition from the super- 
fluid to the normal phase. These features can be seen in 
Fig. QJa) where 5/i = 0.4t and 5/i = 0.7to correspond- 
ing to P w 0.12 and P ss 0.34, respectively. Similar 
spatial modulations have been recently found also in di- 
lute population imbalanced mixtures without an optical 
lattice [H, [H, H3| , however they have not yet been ob- 
served in the current experiments [H, 0, HI 0] ■ -"- n contrast 
to our BdG results, the LDA type methods exclude the 
possibility of order parameter modulations and therefore 
they fail to produce such a spatially modulated super- 
fluid phase which is one of the possible candidates for 
the ground state. 

In the recent theoretical works on dilute population im- 
balanced mixtures without an optical lattice, such spatial 
modulations have been suggested as signatures for the 
FFLO type superfluidity by some authors [H, [l9| and as 
finite size effects by some others [2(| ■ Here we remind 
that the FFLO type superfluidity is characterized by the 
formation of Cooper pairs with nonzero center-of-mass 
momentum, in contrast with the BCS type superfluid- 
ity where Cooper pairs have zero center-of-mass momen- 
tum [22j. Therefore in two- and three-dimensional sys- 
tems it is an open question whether these spatial modula- 
tions are related to the FFLO superfluidity or are simply 
finite size effects. However we also remind that the exact 
ground state phase diagram of one-dimensional systems 
have been recently calculated [H, [H, [2f| [26[ snowing 
that the superfluid phase has FFLO structure in trapped 
as well as infinite systems. 

In Fig. [TJb) we show that the unpaired fermions are 
pushed away from the trap center towards the trap edges 
and they have a maximum at the position where A^ 
changes sign. This is because spatially bound Andreev 
type states form around the nodes of A^, and the oc- 
cupation of these bound states is different for f and J, 
fermions [l8j |. Since /jl^ > /i| when > iVj_, the | 
fermions mostly occupy these states leading to the single 
peak structure. This feature is in good agreement with 
the recent experiments on dilute population imbalanced 



mixtures without an optical lattice [j], @, H, 0] . However 
in contrast with the trapped mixtures without an optical 
lattice, both A^ and pi have C4 symmetry which is con- 
sistent with the underlying symmetry of the square lat- 
tice. Here we notice that the LDA type methods always 
produce results with rotational symmetry and therefore 
they are not strictly applicable to optical lattices. Having 
shown that the ground state phases of low density mix- 
tures in optical lattices are qualitatively similar to those 
of the dilute mixtures without an optical lattice, next we 
discuss the high density mixtures. 
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FIG. 2: (Color online) We show the order parameter Ai (on 
the left, in units of to) and population difference pi = rii^ — 
m t i (on the right, per lattice site) for the high density case on 
a two-dimensional square lattice with 50a x 50a sites. Here 
the chemical potentials are such that (a,b) 5/i = OAto; (c,d) 
5/i = 0.5to; (e,f) S/i = 0.6io; and (g,h) 5/i = 0.7irj- 

(II) High density mixtures: In Figs. [21 and [3] we show 
the superfluid order parameter Aj and the population dif- 
ference per lattice site pi for the high density case where 
N « 1570. When U = and iV T = N t = N/2, the 
maximum filling of this case corresponds to a fully-filled 
band with rii^ = 1 near the trap center. For such high 
densities the ground state phases are very different from 
those of the low density systems as can be seen in Figs. [2] 
and To understand these ground state phases of pop- 
ulation imbalanced mixtures, it is again very illustrative 
to first discuss the population balanced case. 
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For a weakly attracting population balanced mixture 
with Uq = 3to and 8fx = 0, we find that Aj = around 
the trap center for distances |rj| < 4a. This signals the 
band insulator phase characterized by a fully-filled band 
where n^-f = n^i = 1. However since m t -\ = rii t i < 1 
away from the trap center, Aj becomes finite signal- 
ing a transition from the band insulator to the super- 
fluid phase. The maximum Aj occurs around |r;| « 16a 
where m -\ — n%,± = 0.5 corresponding to a half-filled 
band. This is purely a density of states (Dj) effect since 
Aj oc tQe~ 1 /( u ° Di '> and Z3j has a maximum exactly at 
half-filling due to particle-hole symmetry of the Fermi- 
Hubbard Hamiltonian. For longer distances |ri| > 16a 
away from the trap center, A.; gradually decreases until 
it eventually vanishes for |ri| > 22a where n^j = m ^ ~ 0. 
These features can be seen in Fig. 0(a) and they are 
very different from those of the low density case shown 
in Fig. [Ha). In contrast to our BdG results, the LDA 
type methods fail to describe the band insulator region 
with unit filling because the density profiles do not vary 
smoothly as a function of |ri|. 
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FIG. 3: (Color online) We show (a) the order parameter Ai (in 
units of to) and (b) the population difference Pi = n^f — n^j, 
(per lattice site) for the high density case as a function of 
distance x (in units of a) from the trap center. Here y = 0a. 

In the case of population imbalanced mixtures, we 
find that Aj modulates in the radial direction towards 
the trap center for low imbalance as shown in Fig. [2K a ) 
because Aj is a more slowly decreasing function of |ri j 



towards the trap center than towards the trap edges 
when 5 fx = 0. However Aj also modulates towards 
the trap edges with increasing imbalance as shown in 
Figs. EJc) , [21(d) and [21(e). Characteristic features of these 
spatial modulations are similar to those of the low den- 
sity systems and they can be seen in Fig. [31(a) where 
Sfj, = 0.4to, SfJ* — 0.5to, 3fJ- = 0.6to and Sfx = 0.7io cor- 
responding to P w 0.017, P m 0.058, P m 0.090 and 
P rs 0.12, respectively. Therefore high density mixtures 
in trapped optical lattices are also good candidates for 
observation of such exotic superfluid modulations. Fur- 
ther increasing the population imbalance gradually de- 
creases Aj as shown in Fig. [31(a), until it vanishes en- 
tirely beyond a critical imbalance signaling a transition 
from the superfluid to the normal phase. 

In Figs. [U(b) and[3Jb) we show for low imbalanced mix- 
tures that the density of unpaired fermions has a single 
peak at the position where Aj changes sign. However 
since Aj also modulates towards the trap edges for in- 
termediate imbalance, the unpaired fermions have dou- 
ble peaks in their density as shown in Figs. [Hd), |U(f) 
and (3](b). Furthermore since Aj vanishes with further 
increase in imbalance, these two peaks merge leading to 
a single peak which is shown in Figs. [21(h) and [2Kb). No- 
tice that similar to the low density case both Aj and pi 
have C4 symmetry which is consistent with the underly- 
ing symmetry of the square lattice. 

Conclusions: To conclude we used the BdG method 
to analyze the ground state phases of population imbal- 
anced fermion mixtures in harmonically trapped optical 
lattices. First we showed that the phase structure of low 
density mixtures in optical lattices are qualitatively simi- 
lar to those of the dilute mixtures without an optical lat- 
tice. Then we discussed high density mixtures and found 
qualitatively different results. In both cases we found 
that the superfluid order parameter modulates spatially 
but it is an open question whether these modulations are 
related to the FFLO superfluidity or are simply finite 
size effects. Lastly we compared our BdG results with 
the LDA ones and argued that the LDA type methods 
are not sufficient to describe especially the high density 
mixtures in harmonically trapped optical lattices. 
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